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Abstract 

Different aspects of bozonization algorithm proposed in p, || are tested by nu- 
merical simulations of a one dimensional toy model. 

1 Introduction. 

Computer simulations in lattice models with fermions meet serious difficulties due to 
grassmanian nature of fermionic variables. The method mainly used so far was based on 
the Hybrid Monte Carlo algorithm |1|. Considerable progress was achieved in this direc- 
tion but with the present computer facilities it still does not allow efficient calculations 
in models with dynamical fermions. 

An alternative approach was put forward by M. Lusher 0, [3], who proposed to calculate 
a fermion determinant replacing it by an infinite series of boson determinants. However 
at present the efficiency of this algorithm is comparable with the efficiency of Hybrid 
Monter Carlo method mainly due to the large autocorrelation time. (For a recent review 
of different approaches to the problem see M ) . 

So at the moment the problem is still open and alternative approaches to bozonization 
should be investigated. In papers |, | a representation of a D-dimensional fermion 
determinant as a path integral of a (D+l)-dimensional Hermitean bozonic action was 
proposed. This construction provides a new algorithm for numerical simulations of lattice 
QCD with dynamical quarks, although it is generally applicable to any model whose action 
is quadratic in Fermi fields and the matrix of quadratic form is positive definite. To study 
various technical aspects of a simulation procedure we carry out in the present paper a 
detailed study of a one dimensional toy model. In the original paper [f| bozonization 
procedure was based on the constrained bozonic effective action. Later it was observed 
H that using a freedom in the choice of effective bozonic action one can get rid off the 
constraints which seems to be more convenient for Monte-Carlo simulations. In this paper 
we use the unconstrained version. 



2 Effective bozonic action. 

Let us consider a D-dimensional lattice fermionic model with the action: 



S f = a D J2^(x)(B 2 + m 2 )^{x) (1) 

X 

where x numerates the sites of a D-dimensional Euclidean hypercubic lattice, and B 2 is 
some positive bounded operator. (We shall assume also that the operator B is Hermitean, 
although it is not really necessary.) 

We introduce (D+l)-dimensional bozonic fields 4>(x,t) which have the same spinorial 
and internal structure as ip(x). 

The extra coordinate t is defined on the one dimensional chain of the length L with 
the lattice spacing b: 

L = Nb ; < n < N (2) 

We assume that b || B ||<C 1. 

The determinant of the operator B 2 + m 2 can be presented as the following bozonic 
path integral || || : 

Det(B 2 + m 2 ) = [ e'^D^pDip = lim / e - Sb[M D(j)*D(j)Dx*Dx (3) 

where: 

N-l 

S b [4>, X ] = a D Y.b E [-^ 2 (C +1 (x)0„(x) + h.c. - 2<t>* n {x)j> n {x)) + 

x n=0 

+b-\t<f)l +1 (x)B<t) n {x) + h.c) + -{<j)* n+1 (x)B 2 <t) n {x) + h.c) + 

+v / Zexp{-mim}(x*(x)(m + iB)<p n (x) + h.c.) +— — ^2x*( x )x( x ) (4) 

Zm x 

and the free boundary conditions in t are imposed: 

<j) n = ; n < ; n > N (5) 

The bozonic D-dimensional fields xi x ) have the same spinorial and internal structure as 
the fields ip(x). 

The bozonic action (Q) is a linearized version of the expression in the exponent of the 
following integral: 



I = fexp{a D J2b E [b- 2 (€liexp{-iB a b}<p a n + h.c. - 2<CC) - 
-VLexp{-mbn}(x a *(m + zB a )^ + h.c)]-^Y,X a *X a }D4>*D4>Dx*Dx (6) 



where instead of x-representation we used a basis formed by the eigenvectors of the oper- 
ator B, B a being corresponding eigenvalues. 

Indeed, taking into account that b \\ B ||<C 1, we can expand the expression in the 
exponent of the integrand (@) in a Taylor series. Keeping only the terms, nonvanishing in 
the limit b — > 0, we get the expression (||). For a finite b the difference between (|3p and 
(|) is of order 0(b 2 \\ B \\ 2 ). 

By changing variables 



cxp 



{-iB a nb}<p a n ■ C - exp^nfe}^ (7) 



we can rewrite the eq. @ as the gaussian integral over <fr with the quadratic form which 
does not depend on B a . To calculate this integral it is sufficient to find a stationary point 
of the exponent. 

For small b the sum over n can be replaced by the integral and the equations for the 
stationary point acquire the form: 

Q - VZx a (m - iB a )e' {m ^ Ba)t = (8) 

^(0) = 4> a {L) = (9) 

(One can show that replacing the sum by the integral also introduces corrections of order 
0(b 2 || B || 2 .) 

The solution of these equations is: 

m — iB a v L ' 

Substituting these solutions to the integrand ,we get: 

f ( x a *x a 

lim/ = yexp{-E^ T ^ 



In,,/ = fexp{-J2 \ * 2 (1 ~ 2e~ mL cos B a L + e- 2mL )}D X *D X (11) 



Therefore 



lim I = det(B 2 + m 2 ) (12) 

b— >0,L— »oo 



The equality (||]) is proven. For finite h and L this equation has to be corrected by the 

tjGrms" 

0(b 2 || B \\ 2 )+0(e- mL ) (13) 

In the next section we apply this construction to simulations of a one dimensional 
model. 

3 Numerical simulations for free fermions on the 9 1 
lattice. 

Some technical issues of the new algorithm can be observed by performing numerical 
simulations on ID lattice for free fermions. We used the following action: 

S f = J2r k (-d 2 +m 2 )i> k (14) 

k 

where ip*,ip are anticommuting Grassman variables, d is symmetrical lattice derivative 
and m is a mass. The lattice spacing a has been set equal to 1 for convenience. One can 
rewrite the operator in the quadratic form of (0) as follows: 

- d 2 + m 2 = B 2 + m 2 (15) 

where B = {id) is hermitian matrix. So in accordance with the discussion above the 
path integral over Grassmanian variables ip*,ip can be approximated by the path inte- 
gral over bozonic fields 0, \. The corresponding bozonic theory with the multiquadratic 
action 5&[0, x] ( see ec L-(|])) can be simulated straightforwardly using local heatbath and 
overrelaxation algorithms. 



We note that for the model (|14]) the action S b [<f), x] can be rewritten in the form: 

S b [<f ) ,x} = S[p,(3} + S[a, 7 } (16) 

where p, j3 and a, 7 are real and imaginary parts of the fields <ft, \: 

cf) = p + 10 

X = f3 + t7 (17) 

Hence the path integral over the fields 0, \ allows the following factorization: 

Z b = f e~ Sb[ ^ x] D(j)*D(t)Dx*Dx = [ f e-~ s[pA D pD (3~\ 2 (18) 

It makes possible to simulate the theory with the real fields p, (5 and the action S[p,/3], 
accelerating the calculations by the factor of 2. The resulting partition function must be 
squared. 

Let us firstly consider the updating of the field p. If all field variables except p at 
point (x,t) are kept fixed (where t numerates points in auxiliary dimension), the action 
assumes the following form: 

S b [p(x, t)] = A(p(x, t) - E(x, t)f + const (19) 

where A is a positive constant and E(x,t) is an easily calculable vector. A local update 



p(x, t) -► p(x, t) = uE(x, t) + (1 - u)p(x, t) + ^ A — V (20) 

where r] is a gaussian random number of unit variance, fulfills detailed balance for any 
0<lu< 2 0. 

In our test we used hybrid overrelaxation algorithm, which consists in the mixing of 
heatbath (u = 1) and overrelaxation (u = 2) sweeps with a ratio 1 : N or ||. It is believed 
that this algorithm has a dynamical critical exponent z m 1 if N or is proportional to the 
correlation length £. 

Generally we subdivided the lattice into 2 sublattices, coloring each site according to 
the function 

c( x ,t) = (-iy 

One can see that E(x, t) does not depend on the field variables at the same t and sites of 
the same colour do not interact with each other. We updated one colour after the other. 
Analogously, if all field variables except (3 at point x are fixed, the action assumes the 
form: 

S b [p(x)] = C([3(x) - D(x)f + const (21) 

where C is a positive constant. Since field j3 is interacting with p(x,t) for all t, the 
computation of the vector D(x) requires an effort proportional to N, i.e. a j3 field update 
is almost as expensive as the updating of p field. 

In our implementation, the iteration is made up of one p heatbath sweep, Ng r overre- 
laxation sweeps, one j3 heatbath sweep and N&. overrelaxation sweeps. After each iteration 
the following function was measured: 

K=< y £Wk> (22) 

k 



To measure the function ( p2|) using the bozonic approximation of the fermionic theory 
with the action (|14D, one must rewrite (|22|) as a correlation function of p,/3 fields with 
the measure defined by S[p, /?]. Let us denote: 

Z/ = / 'e~ s fDilj*Dtjj (23) 

Z h = J e-^DpDfi (24) 

In the section 2 it was proved that 

Z/ = Z 2 + 0{b 2 || B || 2 ) + 0{e~ mL ) (25) 

Then we can derive: 

n- l 8 z- l 5 z-- l 8 z- 1 < ^M1 > (26) 



and 






We used the expression fl27|) as a bozonic approximation to the function (|22|). One can 
see that the accuracy of this approximation is also given by the expression (|13|) . 

In table 1 the autocorrelation time dependence is displayed for several updating 
schemes. Autocorrelation times were measured for the function (|27|) using the method 
proposed by Socal 0, namely 

*-(*•> = 5 + !§§ (»> 

with 

i n—i 

C[i) = : £ {K k - TZ) (K k+l - H) (29) 

where the n is a number of iterations and M chosen so that r int <C M <n. An estimate 
for the error of Ti nt is given by 



2(2M+1) 

a r snt = r int (30) 



One can see that overrelaxing (3 field does not decrease autocorrelation time substan- 
tially (it even may increase T int in units of CPU time). Contrary, overrelaxing p field 
improves the autocorrelation behavior. When adding more overrelaxation sweeps for the 
given sets of parameters, r int in CPU units starts to rise again. 



Updating 


n b 


n n t{T^b) 


Hh 


0.526(10) 


56(11) 


HOh 


0.534(7) 


20(3) 


HhOo 


0.544(9) 


22(4) 


HOhO 


0.546(6) 


10(1) 


HOhOo 


0.546(6) 


9(1) 


HOhOO 


0.544(4) 


4.3(3) 


HOOhOO 


0.543(3) 


2.8(2) 


HOOOhOO 


0.541(3) 


2.2(2) 


HOOhOOOo 


0.541(2) 


1.6(1) 



Table 1: Autocorrelation times in [1/iteration] units on 9 1 lattice for m = 4, N = 100 and 
b = 0.015. The letters in the first column give the type and order of sweeps used per iteration, 
where H is a p heatbath, O is a p overrelaxation and h and o are the (3 updates. The exact value 
of n is 0.5457. 



Now we discuss the systematic error and autocorrelation behavior of the algorithm. 
We measure the function H b (see (|27|)) and r int {lZb) for different sets of parameters b and 
N. As it is demonstrated in section 2, the systematic error is given by: 



where: 



A = A! + A 2 , 

A 1 = Cb 2 + Fb 3 + 0{b 4 ), 
A 2 = De- mbN , 



(31) 

(32) 
(33) 



and C, F, D are some constants. 

To control the systematic error one can fix the parameter mbN and perform calcu- 
lations for different values of b. We choose m = 4 , mbN = 8 . The results are shown 
on Fig.l where the function lZb{mb) ( ^7|) is plotted (the statistical errors are small). The 
horizontal line corresponds to the theoretical value (|22|) . It is seen that the results con- 
verge to the theoretical value as mb decreases. The fit of expression (|32|) for Ai gives 
C ~ 17 and F ~ —43. The systematic error 0(b 2 ) is partially compensated by the 0(b 3 ) 
error. 

On Fig. 2 the dependence of the autocorrelation time for TZ^ against mb is plotted. In 
these and all further measurements HOhOO scheme is used. From these data we get that 
Ti n t(TZb) ~ 0.25/mb when mb — > 0. 

Now let us fix the parameter b = 0.01 making Ai less than 0.002 and measure TZf, 
and T int (TZb) for the different number of points in auxiliary dimension N. On Fig. 3 the 
dependence of 7^ against A" is shown. The horizontal line denotes the theoretical value 
of (F2"2"|). From these data we get D < 2.5. 

On Fig.4 the dependence of T int (7Zb) against A^ is plotted. The autocorrelation time 
grows quadratically when A^ increases, but the proportionality factor is very small: 
~ 10~ 4 . For large A^ the autocorrelation time can be decreased by adding more overrelax- 
ation sweeps of the p field. To investigate the autocorrelation behavior of the algorithm for 
the cases of practical importance one needs to study the models of larger dimensionality 
with the gauge fields involved. 

In conclusion we investigate the slowing down of the algorithm at small values of m. 
We fix b = 0.015 and mbN = 6 making the systematical error constant. The results are 
shown on Fig.5, the fit of these data gives r int (TZ b ) w -^, C w 3.3, a ~ 2. 
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4 Discussion. 

We performed the first simulations for Slavnov's algorithm on ID lattice for free fermions. 
It was shown that correct and accurate results can be obtained with a reasonable size of 
lattice in auxiliary dimension. We are going to extend this simulations to larger lattices 
taking into account interaction with the gauge fields (in progress). 
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